function [A, B, Phi] = WAIA_step1(I,Delta,W)
    [I_row, I_col] = size(I{1,1});
    A = nan(I_row, I_col);B = A;Phi = A;
    S = AIA_step1_S(Delta,W);
    for r = 1:I_row
        for c = 1:I_col
            R = AIA_step1_R(I,Delta, r, c,W);
            R_nan = isnan(R);
            if sum(R_nan) == 0
                X = S \ R;
                A(r,c) = X(1);
                Phi(r,c) = atan2(-X(3), X(2));
                B(r,c) = sqrt(X(2)^2+X(3)^2);
            end
        end
    end
end

function S = AIA_step1_S(Delta, W)
%     F = length(Delta);
    S11 = sum(W(:));
    S12 = W.*cos(Delta);S12 = sum(S12(:));
    S13 = W.*sin(Delta);S13 = sum(S13(:));
    S22 = W.*cos(Delta).*cos(Delta);S22 = sum(S22(:));
    S23 = W.*cos(Delta).*sin(Delta);S23 = sum(S23(:));
    S33 = W.*sin(Delta).*sin(Delta);S33 = sum(S33(:));
    S = [S11, S12, S13;
         S12, S22, S23;
         S13, S23, S33];
end

function R = AIA_step1_R(I,Delta,r,c,W)
    F = length(I);
    R11 = 0;R12 = 0;R13 = 0;
    for f = 1:F
        delta = Delta(f);
        i = W(f).*I{f,1}(r,c);
        R11 = i + R11;
        R12 = i*cos(delta) + R12;
        R13 = i*sin(delta) + R13;
    end
    R = [R11, R12, R13];
    R = R';
end